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The dynamic behavior of microtubules in solution can be strongly modified by interactions with 
walls or other structures. We examine here a microtubule growth model where the increase in size 
of the plus-end is perturbed by collisions with other microtubules. We show that such a simple 
mechanism of constrained growth can induce ordered structures and patterns from an initially 
isotropic and homogeneous suspension. First, microtubules self-organize locally in randomly oriented 
domains that grow and compete with each other. By imposing even a weak orientation bias, external 
forces like gravity or cellular boundaries may bias the domain distribution eventually leading to a 
macroscopic sample orientation. 



I. INTRODUCTION 

Biological processes like cell division, transport of cer- 
tain organelles, morphogenesis and organization in the 
cell are mediated by rod like structures known as micro- 
tubules, which form various arrays, radial spindles, par- 
allel and antiparallel bundles 0, 0, . The microtubule 
self-assembly in living organisms is regulated by different 
factors: microtubule-associated proteins (MAPs) which 
stabilize, destabilize and crosslink microtubules 0,01 di- 
verse kinesin-like motor proteins, which organize and link 
microtubules, 7-tubulin ring complex which serves as a 
template for nuclcation sites for microtubule polymer- 
ization in centrosomes |(J • These factors combine with 
the physical and chemical properties of the solution to 
determine, by mechanisms not yet well understood, the 
spatial structure and the orientation of the microtubules. 

Each individual microtubule is a highly dynamic self- 
assembled rod, which is permanently growing or shrink- 
ing. This ability for being in an everlasting state 
of changing length as won microtubules the name of 
"searching devices" for specific targets in the cell 0,0. 
A key property allowing for this bistable state is the dy- 
namic instability Q. Due to conformational asymmetry 
of the constituting microtubule subunit, the heterodimcr 
a,/3-tubulin, a microtubule has a polar structure and it 
grows at a plus-end and shrinks at its mmns-end. The 
speed of growth at the plus-end is not constant, but 
rather intermittent. The elongation of the plus-end is 
stochastically alternated by the abrupt shrinking, in a 
process of unidimensional diffusion |J. Such dynamic 
behavior is attributed to the complex, two-stage assem- 
bly of the plus-end, implying the internal hydrolysis of 
the GTP in a tubulin dimer. Within the cap model 
0, HI, , the tubulins added to the growing plus-end are 
not hydrolyzed, thus having configurations favorable for 
the microtubule assembly. They form a cap protecting 
the plus-end from the disassembly and shrinking. Tubu- 
lins incorporated into microtubules are capable for hy- 
drolysis. During conversion of the GTP of /3-subunit to 
the GDP, the tubulin heterodimer undergoes the confor- 
mational change that destabilize the tubular structure 
of a microtubule and favors shrinking The rate of 



hydrolysis compete with the polymerization rate of the 
plus-end. Some factors like local fluctuations of tubulin 
concentration near the plus-end may perturb the balance 
thus leading to the loss of the protecting cap and provok- 
ing rapid shrinking. Since the density fluctuations have 
a stochastic nature, the growth to shrinking transitions 
(catastrophes) and the inverse resume of growth (rescue) 
are also statistically distributed. 

The catastrophes can also be induced by factors other 
than concentration fluctuations. During the growing 
state, microtubules interact with each other as well as 
with different cellular structures. The mechanical force 
opposing the growth of the plus-end may induce catas- 
trophes |l2l Experiments involving growing mi- 
crotubules and different immobile obstacles and barriers 
have shown that the opposing force reduces the growth 
velocity 0,0]. Thus, the dynamic instability coupled 
with the mechanical force may influence the microtubule 
length and thus induce a spatial organization. Indeed, 
the catastrophe rate is expected to be higher near cel- 
lular boundaries and lower in the cytoplasm |l2j . The 
boundaries may also induce the orientation preference, 
since the catastrophe rate in the direction perpendicu- 
lar to the boundary would be higher then that along the 
boundary. Recent in-vivo work suggests that the micro- 
tubule dynamic instability is altered during preprophase 
band formation |l(j| . Microtubule reorientation is accom- 
panied by the increase of the catastrophe frequency and 
growth rate, while the rescue frequency and shrinkage 
rate remain unchanged. The gradients always present in 
living cells can also play a role of the "effective" bound- 
aries and provoke the microtubules ordering |17| . In par- 
ticular, the gradients of the energy dissipation, concen- 
tration of tubulin and associated proteins may result in 
spatial anisotropy of the growing and shrinking speeds 
which can lead to a self-organization. 

Another important observation concerns the inter- 
microtubules interactions [l8| . Encounters between cor- 
tical microtubules affect their dynamic behavior. Steep 
contact angles of microtubules collisions provoke catas- 
trophes more often than the shallow contact angles, while 
the microtubules with close angles have shown a tendency 
to zippering. 

These observations suggest that the spatial self- 
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organization of microtubules might be induced by the 
coupling of the catastrophe events and the configuration 
of neighboring microtubules. Different examples of spon- 
taneous self-organization of dynamic microtubules have 
been reported in the literature yjj, l20l l2ll l'22l |23| . One 
of them |23, describes the formation of spatial struc- 
tures in the in-vitro solution of microtubules that start 
growing from seeds distributed homogeneously. The re- 
sulting pattern strongly depends on the direction of the 
gravitational field, where microtubules are organized as 
highly aligned strips. Interestingly, the observed struc- 
tures do not appear in weightlessness conditions. The 
authors conclude from their observations that the Earth 
gravity only triggers the symmetry breaking and does 
not affect individual microtubules. The concentration 
of microtubules even coupled with gravity is not suffi- 
cient to provoke the system orientation due to excluded 
volume effects as in usual liquid crystals p3. |25| : the 
shaking or mixing of the sample irrevocably destroys the 
pattern. Since the in-vitro preparation does not contain 
any molecular motors or MAPs, the self-organization in 
stripes is attributed to the dynamic nature of the mi- 
crotubules: the pattern formation disappears when the 
dynamic instability is inhibited by the addition of taxol. 
Similar patterns are formed under magnetic fields |2lj |. 
The in-vivo observations of the reorientation of cortical 
microtubules in parallel arrays [2^, I23L l2rj [2?J are also 
ascribed to dynamic instability: the depolymerization of 
disordered microtubules is followed by the repolymeriza- 
tion into an ordered array. However, the in-vivo self- 
organization can be regulated by MAPs or motor pro- 
teins 0, or it may be a result of the simultaneous 
A microtubule is modeled as a rigid, oriented rod 
which shrinks at its minus-end and grows at its plus- 
end. Instead of dealing with a fluctuating rate of growth 
and shrinking, we rather consider smooth, coarse-grained 
properties, namely the average speeds of growth and 
shrinking. Without obstruction, in a free environment, 
the plus-end of a microtubule grows at constant speed i> + , 
while its minus-end shrinks at constant speed v~ . When 
the plus-end encounters another rod, it stops, but the rod 
continues to shrink at its minus-end, with speed v~, and 
as a result, the overall length of the rod decreases. The 
rod resumes its growth as soon as its plus-end is no longer 
blocked by its neighbor. Altogether, the plus-end expe- 
riences an environment dependent intermittent growth, 
and the minus-end a constant motion at speed v~. The 
rod disappears if its length decreases to zero during the 
shrinkage phase. The total number of rods is not fixed, 
but is maintained by a permanent injection rate of new 
rods at random positions, with random orientation and 
zero length. 

We found that this mechanism of the constrained 
growth can lead itself to the spontaneous alignment of 
microtubules in an initially isotropic and homogeneous 
array. The emergence of a local orientational order is 
reminiscent from a natural selection process. Whenever 
some local anisotropy builds up, the survival rate of the 



action of many factors. 

In this paper we explore the consequences of the sim- 
plest possible physical hypothesis for explaining micro- 
tubule orientation from a coupling mechanism between 
growth and inter-tubule interactions. In the next section 
we set the foundations for the physical model based on 
the observation that the inter-microtubules collisions in- 
crease the rate of catastrophes [l8| ■ In section lnTl we show 
numerically that this mechanism alone leads to the ori- 
entation of microtubules in aligned stripes. A theoretical 
discussion of our main results is presented in section llVl 
and our findings are summarized in the conclusion. 



II. THE MODEL 
A. A kinetically constrained growth model 

We propose a model based on the assumption that 
the assembly dynamics of a particular microtubule is 
influenced by the configuration of the surrounding mi- 
crotubules. Our purpose is to exhibit the simplest pos- 
sible model of constrained growth inducing some self- 
organization of the microtubules. Thus, in our model, we 
neglect all the other possible physical mechanisms that 
could play a role in the co-alignment, such as collision 
induced turnover and dislocation |22j |. as well as the ex- 
cluded volume interactions leading to ordering in usual 
liquid crystals. The incorporation of these factors would 
facilitate and enhance the alignment, 
neighboring rods changes and becomes orientation de- 
pendent. The rods which from the start have picked up 
the dominant orientation are likely to outlive rods with a 
different orientation, and this creates conditions favoring 
the population of rods with the "correct" orientation at 
the expense of rods with "incorrect" orientation. In our 
model, the rods cannot change their orientation but the 
permanent injection of young, randomly oriented rods, 
leaves to the system the possibility to reorganize and to 
tune up to a change of external conditions. 

We implemented numerically this constrained growth 
model. The simulations clearly show a trend towards 
some local organization and ordering, with the formation 
of well defined anisotropic domains. 



B. Numerical implementation of the constrained 
growth model 

In our numerical implementation, in two dimensions, 
each rod is characterized by its length, its position, and 
its orientation. The orientation (one angle) and the po- 
sition (two coordinates) are set when the rod is injected, 
and do not change until the rod eventually disappears. 
The length of each rod evolves with time, starting from 
zero shortly after the injection. All rods are packed in 



FIG. 1: Snapshots of the system of rods at various successive times (time increases from a to d). Growing rods are drawn in 
blue, and shrinking rods in red. 



a square box of side L s , subject to periodic boundary 
conditions. 

Depending on their dynamic state, the rods belong to 
two categories: shrinking rods (s-rods) are the rods that 
the local environment prevents from growing further (ki- 
netic constraint), while growing rods (g-rods) are free to 
grow. 

The rule for the time evolution of the rods is the non- 
crossing displacement. Updates are done every time in- 
terval At. The minus tip of each rod is shortened by the 
amount v~At, while the plus tip attempts a move for- 
ward by v + At. If this move can be done without crossing 
any other rod, the move is accepted and the rod grows. If 
the rod was in a s-rod state, it converts to a g-rod state. 
Otherwise, the move is rejected, and the rod switches to, 
or stays in a blocked s-rod state. As a result, the length 
of a rod after each step, either increases by (v + — v~)At, 
or decreases by At. 

In the present case, we chose the values L s = 100 and 
At = 1. The speed v + ranges from 0.5 to 2.1 and v~ = 
0.3. In the discussion, we make a frequent use of the 
speed ratio a: 



defined as the ratio between the speed of shrinkage in the 



s-state v~ and the speed of growth in the g-state v + — v ~ . 
The corresponding values of a used in the simulation lie 
in the interval 0.17 to 1.5. The value for the injection 
rate Qi is about 1000 new rods per unit of time. There 
are typically a few thousand rods (5000 to 6000) at a 
time in the simulation box. 



III. MAIN RESULTS 
A. Spatial organization: domain structure 

The simulation starts with a set of rods of zero length. 
The kinetic constraint concerns a vanishingly small num- 
ber of rods at the early stage of the system evolution. 
As the number and the length of the rods increase, the 
amount of packing gets larger, and the kinetic constraint 
forces a significant fraction of rods into a blocked, shrink- 
ing state. This transient regime recedes to a quasi- 
stationary regime in which the ratio of s-rods and g-rods 
seems to remain approximately constant. 

Then, the numerical simulations shows clearly the slow 
emergence of oriented domains, or bundles, of nearly 
aligned rods. These domains show very rough, ragged 
and sharp boundaries, much as crystallites. Thus, they 
look quite different from the domains arising in the usual 
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FIG. 2: The dominant angles O, 0; as a function of the time, 
as found by minimizing equations and |J3 



phase transitions and coarsening situations, where a fi- 
nite bending elasticity creates a smooth variation of the 
order parameter, at the vicinity of a domain wall. The 
domains seems to be randomly oriented, and the isotropy 
of the system is recovered only on length scales larger 
than the size of the domains. 

In the final stages of the simulation, the average do- 
main size L(t) is still a slowly growing function of the 
time, and eventually becomes of the same order of mag- 
nitude as the size of the system L s . This prevents reach- 
ing an asymptotic finite value of L(t), associated to a 
truly stationary distribution of the rod lengths and ori- 
entations. Well known examples of such coarsening dy- 
namics are characterized, for instance, by a power-law or 
a logarithmic behavior of L with t [2^| . In the latter case, 
the possibility to discuss the system properties in term 
of quasi-stationary solutions remains. 

Figure illustrates how the rods self-organize with 
time. At first, the population of rods is isotropic, ex- 
cept for small fluctuations inherent to the random ini- 
tial position and orientation distribution (Figure QJi). 
These pre-existing heterogeneities grow into small bun- 
dles, whose distribution still remain seemingly isotropic 
on large scales (Figure^). Then, larger bundles emerge 
at the expense of many other smaller bundles, bound 
to disappear (Figure^). Finally, the typical size of the 
larger bundles becomes comparable to the size of the sim- 
ulation box (Figure QJi) . The absence of bending mod- 
ulus, and the presence of ragged boundaries, forbids a 
mechanism based upon domain walls motion. Such a 
competitive growth of the domains is slow. 

To quantify the degree of local ordering, or "polariza- 
tion", of the system, we introduce a dominant angle 0, 
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FIG. 3: Anisotropy ratios Si and S vs time t. 



which maximizes a cost function a: 



a(9) 



— V cos 2 (fi< - 9) 



cos 2 (n-0), 



(2) 



where Qi is the orientation (angle) of the i rod, and 
the sum runs over nt, the total number of rods present in 
the system. Thus, a is defined as the "ensemble average" 
over the population of rods at time t, denoted with an 
overline . . . We call order parameter the value of the 
maximum s = cr(0). The cost function can be expanded 

as a (9) = cos 2 Q cos 2 # + sin 2 fl sin 2 9 + 2 sin Q cos fl sin 2 9. 
Then, the stationarity condition da{9) / 89\q = leads to: 



tan 29 = 



i2fl 



cos 2fi 



(3) 



This equation has always four solutions, two correspond- 
ing to the maxima O max and max + tt 5 and the other 
two, to the minima m i n or m i n + tt, and Q m ax and 
©min are mutually orthogonal. 

A scaled anisotropy parameter S may be defined as: 



^(Qmax) ~ ^(Qmin) 

c(S max ) + cr(e min ) ' 



(4) 



The quantity S is for a population of isotropically ori- 
ented rods, while it is 1 for a population of perfectly 
aligned rods. The parameter S makes it possible to quan- 
titatively assess the amount of ordering in the system. 
Because both parameters s and 5 turn out to fluctuate 
strongly with time, we introduce a more stable parame- 
ter, where each rod i contributes according to its length 
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FIG. 4: The anisotropy ratio Si vs time, for three increasing 
biases e in the angular distribution of the rods orientation. 



FIG. 5: Different rod histories, showing the length l(t) func- 
tion of the time t. The curves are similar to random walks 
with an absorbing boundary condition at I = 0. 



where the long rods participate more than the short 
ones. The dominant angle associated with this parameter 
obeys: 



tan 29, = 



I 2 sin 2n 
I 2 cos 2fl ' 



(6) 



and the anisotropy ratio Si defined as in Eq. 0] with a 
replaced by ct,. 

Figure [3 shows the variation of the dominant angles 
O and 0/ with time, for one set of parameters. At the 
beginning, the system is homogeneous and the distribu- 
tion of angles is isotropic, resulting in a singular and 
noisy function of time. As the system evolves, the or- 
dered structures appear and the angles stabilize around 
their preferred value. It is noteworthy that the 0; curve 
is smoother than the curve, due the stabilizing contri- 
bution of the longest and most stable rods. The plateau 
value is related to the orientation of the dominant bundle, 
and fluctuates from sample to sample. 

The same conclusion can be drawn from the plot of 
the anisotropy ratios S and Si , function of time in Fig- 
ure El Although evolving on the same time scale as S, 
the quantity Si reaches a value closer to 1. The differ- 
ences between the two curves is most certainly due to the 
contribution of the many young, short rods, upon which 
the kinetic constraint has not been acting long enough to 
force them into the dominant orientation. 

The anisotropy ratio and the dominant angles aim at 
quantifying the degree of local ordering in a suspension of 
rod like objects, irrespective of the underlying alignment 
mechanism. In that respect, they make possible a direct 
comparison with other alternative models of microtubule 
orientation. From these curves, one can infer a charac- 
teristic time t* for the emergence of a global orientation 
in the sample, such as, for instance, Si(t*) = 1/2. In 
Figure |3 this ordering time is about a few hundred steps 
(f* ~ 300). 



B. Sensitivity to external stresses 



The constant renewal of the rods, along with the 
growth of the competing domains, confers to the system 
the ability to respond to external perturbations. One of 
our main motivation is to evaluate the sensitivity of the 
ordering to the presence of an external gravitational, or 
magnetic field. Quite similarly, the presence of a hard 
wall is expected to align the nearby domains along its 
direction. 

In order to probe the ability of the rods suspension to 
cope with external constraints, and to monitor its sus- 
ceptibility to a small symmetry breaking, we performed 
several simulations with a slightly biased distribution in 
the orientation of the newly injected rods. Instead of be- 
ing isotropic, we added a fraction e of rods in excess, with 
an angle f2 belonging to a small interval Qo ± 1°. As a 
result, a value e = 0.5% brings about an acceleration of 
the ordering time t* by a factor 1.5, and a value e = 5% 
triggers a three times faster growth of the domains (Fig- 
ure 0J. In both cases, the preferred orientation is clearly 
related to the orientation of the bias £Iq. 

Boundaries and impurities can also induce the align- 
ment. When one of the periodic boundary conditions 
is replaced by a hard wall, the alignment of the rods is 
much faster, and the wall orientation propagates into the 
bulk of the suspension. An identical behavior is observed 
when a rod with fixed length, position and orientation, 
is forced into the simulation box. The rods orient them- 
selves parallel to the guiding rod, and longer guiding rods 
provoke faster ordering. 
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FIG. 6: Scatter plot of the times of growth T and T , as a ' 
function of the total age T — T + + T~ , for a given population 

at time t. Inset: a close-up look at the region of young rods. FIG- 8: Scatter plot of the length (vertical axis) and the age 

(horizontal axis). 
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FIG. 7: Histogram of the distribution of the ages T of a pop- 
ulation of rods, in logarithmic coordinates. The initial decay 
rate is close to 1/T, followed by an exponential decay. 



C. Kinetics of individual rods 

The numerical simulation makes it possible to track a 
single rod as it evolves with time. We observe that during 
its life cycle, a rod can experience many alternating peri- 
ods of growth and shrinkage. A typical microtubule life 
history plot is shown in Figure 03 following a saw-teeth 
curve. 

We call "age" T, the time interval elapsed since the rod 
was injected with zero length in the system. The age is 
the sum of the growing time T + and the shrinking time 
T~ , and the length of the rods can be expressed with the 



help of T+, T~ as 



I = (v+ -v~)T + -v~T~ 
T = T+ + T~. 



or equivalently, 



(v+ -v~)T-l 



T 



v~T + l 



(7) 



(8) 



A typical distribution of both T + and T~ , as a func- 
tion of the age T, is shown in Figure El for a population 
of rods at a given time t (scatter plot), while the inset 
of Figure [B] is an enlargement of this plot in the small 
T region. The values T + and T~ of old rods (large T), 
concentrate near two boundaries, which correspond to a 
length I = in Eq. |SJ For these old rods, which have 
survived many collisions, the growth and shrinkage peri- 
ods compensate almost exactly, and in Eq.0 the length I 
results from the difference between two large quantities. 
By contrast, the young rods (small T) show all possible 
combinations of T + and T~ (inset in Figure EJ- This in- 
dicates that the young rods have not yet been influenced 
by their surrounding. The maximal possible length of a 
rod occurs in the extreme case T~ — and T + = T, thus 
corresponding to a length Z max = (v + — v~)T (the dashed 
line in the inset of Figure [SJ. 

Young rods enjoy a fast growth rate, but many are also 
eliminated quickly. Older rods show a smaller average 
growth rate, but their survival rate increases with their 
age. This is clear from the histogram of the ages, which 
is clearly not exponential, but rather well approximated 
by a power law p s (T) ~ T _1 (Figure [TJ. The relative 
disappearance rate of rods aged T, is p~ 1 dp s /dt ~ T _1 , 
instead of remaining constant, as in the exponential case 
(e.g like for instance the decay of radioactive elements). 
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FIG. 9: Scatter plot of the orientation (horizontal axis) and 
the age (vertical axis) of the rods at three different stages of 
the evolution of the system. The time t increases from top to 
bottom. 



In a sense, the young rods shows "plastic" properties, 
and account for the adaptability properties of the rods 
suspension. By contrast, older rods are expected to show 
more rigidity and persistence from the past history of the 
suspension. We did not find any simple justification for 
the exponent — 1 , which is a numerical finding. This con- 
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FIG. 10: Semilogarithmic plot of the histogram of the num- 
ber of rods, for three different stages of the evolution of the 
system. Straight lines correspond to an exponential decay. 



tradicts a naive argument based on the usual random 
walks, which would predict a T -3 / 2 behavior associated 
to the first return to the origin time distribution. We still 
believe in the analogy, but we attribute this discrepancy 
to the existence of strong correlations between orienta- 
tion, age and life time of the rods. Finally, the age of the 
very old rods is distributed exponentially (Figure 0). 

A scatter plot of the lengths I versus the ages T of 
a population of rods does not support the presence of 
strong correlations between these two parameters (Fig- 
ure HJ. 



D. Distribution of lengths and orientations 

The ordered structures (bundles or domains) effec- 
tively select the rods, keeping only those with an orien- 
tation compatible with the dominant orientation of the 
bundles. The correctly oriented rods collide less often 
with their neighbors than the rods with transverse orien- 
tations, and their "fitness", or survival ability is greater. 

Figure shows three typical scatter plots of the ages 
as a function of the angles. The age of the system in- 
creases from the top to the bottom plot. The presence 
of anisotropic domains manifests itself as sharp peaks 
around a few well defined angles. On the example shown, 
one can see two different domains, respectively around 
25° and 160°. The peaks are duplicated (mirrored) be- 
cause the bundles contain a mixture of two antiparallel 
populations, separated by 180°. 

A typical distribution of lengths regardless to the rods 
orientation is shown in Figure 1101 The distribution is 
exponential except for a very small region of tiny lengths. 
Most of the rods with tiny lengths are very young rods 
injected in the system a few steps ago. They did not 
have time to experience collisions and their distribution 
did not acquire the same characteristics as the old rods. 
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The system not only adjusts the orientation and the 
length of the rods, but it also regulates the total number 
of rods rit- An increase in the rods density leads to a 
higher collision rate among rods, and a higher elimination 
rate. In the early stage of the simulation, in a sparse 
system, the newly injected rods do not meet any obstacles 
and the total number of rods increases sharply. After a 
transient regime, n t evolves very slowly, although it is 
not strictly constant. As a matter of fact, n t slightly 
decreases with the spreading of the dominant orientation 
and the emergence of domains fFigure ITT|l . 



E. The three dimensional case 

Finally, we performed a limited number of simulations 
in three dimensions, in order to check whether the kinetic 
ordering was a special feature of the two dimensional sys- 
tems, or whether it was a generic feature also in three 
dimensions. In this case, our simulation box is a cube of 
size 100 x 100 x 100. In addition to the injection rate 
Qi, to the speeds v + and v~, there is another relevant 
parameter: the diameter, or thickness, of the rods d. 

It turns out that the behavior of the system in three 
dimensions is quite similar to the one observed in two 
dimensions. The initially homogenous solution becomes 
gradually structured into bundles and domains. How- 
ever, the three dimensional system differs by the absence 
of sharp boundaries between domains. The competition 
between the different orientations is not so drastic, since 
bundles with different orientation can interpenetrate if 
the rod thickness is small. Domain walls are more diffi- 
cult to identify, but the main result, i.e. local ordering, 
holds also in three dimensions. 



IV. THEORETICAL DISCUSSION 

We discuss in this part some observed features of our 
numerical simulations: exponential tails in the length dis- 
tribution, collision rates, anisotropy. For this purpose, we 
propose an elementary kinetic theory, and its predictions 
are compared with the numerical simulations. 



A. Ensemble and time averages 

Much like in the usual statistical mechanics, one intro- 
duces two kinds of averages. The time average, denoted 
by brackets (...) corresponds to the mean value obtained 
by the repeated observation of single rods evolving with 
time. For instance, (T + ) is the average growth time of 
a rod. Time averages are accessible through numerical 
simulations. 

The ensemble average amounts to considering the 
whole population of rods at a given time. This corre- 
sponds to an instantaneous "snapshot" of the population 
of rods, and the corresponding average is denoted with 
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FIG. 11: Evolution of the number of rods n t with time, 
and repartition between shrinking rods (n a ) and growing rods 
(n a ). 



an overline . . . : for instance, I is the average length of 
the rods. In practice, the ensemble average consists in 
summing over all the rods, and then dividing by their to- 
tal number n*. The ensemble averages can be computed 
from the simulations, but are also the natural outputs of 
the kinetic theory sketched below. 

The connection between time and ensemble average is 
by no means obvious. In the case of a true stationary 
situation, both averages are expected to coincide. Our 
situation, however is not a full stationary situation, as 
we witness the emergence of unbounded large bundles. 
A slow coarsening dynamics, however, may still exhibit a 
satisfactory agreement between the two kinds of averages. 

Ensemble averages are conveniently handled by means 
of distribution functions. Denoting the length, orienta- 
tion and position of the rods respectively by Z, f2 and r, 
we define Ct(l, fl, r, t) = c g (l, f2, r, t) + c s (l, ft, r, t), where 
c g , c 8 and Ct stand respectively for the distribution of 
the population of growing rods, shrinking rods and to- 
tal number of rods, per unit of surface, at time t, with 
< I < 00, < fl < 2ir and position r. 

Successive integrations over the variables r, I or fi, give 
rise to a hierarchy of distribution functions. In particular, 
the numbers of rods n s and n g are given by: 



n g (t) = drdldttc g {l,n,r,t); 



n s (t) 



drdl dQ,c s (l,Q.,f,t). 



(9) 
(10) 



In what follows, we use a loose notation for the partial 
distribution functions, where the variables which do not 
explicitly appear in c g have been implicitly integrated 
over, e.g. c g {l, Q)dl dQ stands for the fraction of g-rods, 
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FIG. 12: Test of the relations Eq. [Hand Eq. [Bll Circles: 
inverse of the average shrinkage interval (r ~ ) , last term of Eq. 
1181 squares: r.h.s of Eq. 1141 triangles: computed value of p ag - 
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with length between I and I + dl, angle between Q and 
fl + dfi, but located at any position r of the system. 

B. The distribution of lengths 

Following the lines of Appendix EI we restrict our- 
selves to the homogeneous, r-independent case, and as- 
sume that the distributions c g and c s comply with the 
following master equation: 

/ J)t C ff + V 'S1^9 = PsgCs — PgsCg] 

The properties of the distributions mostly depend on 
the injection rate Qi, the speed v = v + — v~ , and the 
speed ratio a — v~/(v + — v~). We introduce the in- 
terconversion rates p sg (£l) and p gs (Q,) between s and in- 
states, and we take care of a possible dependence on the 
direction fi. We find that, for an homogeneous and sta- 
tionary system, the lengths are exponentially distributed 
and verify: 

c = *expH/I(fi)]; (12) 
v 

c s = — exp[-//7(f2)l. (13) 
av 

Such an exponential distribution can be seen in Fig- 
ure El This model accounts for the possibility of 
anisotropic length distributions, via the angle dependent 
function It is possible to find an isotropic, self- 

consistent solution for 7(£!) = 7, but we found also evi- 



dence for an anisotropic self-consistent solution, with a 
non trivial function 7(£1). 

Our predictions for the isotropic solution include a de- 
termination of the rate p sg : 

Psg = (14) 

a determination of the rate p gs : 

= « + | (ir), (15) 

and a self-consistent determination of the average length 
7, as a function of Qi, a, the total number of rods nt and 
the surface S: 

I= V^iyl"(£) 2 - (16) 

The predictions of Eq. El are shown in Figure El an d 
discussed also in Appendix [Bl Agreement is poor for 
short times and it improves for long times. 

The predictions of Eq. 1151 and Eq. 1161 are summarized 
in Figure El The agreement is good for p gs and quali- 
tative for 7 at short times. The prediction becomes poor 
for times larger than t* , associated to the emergence of 
the domains. We believe that the disagreement is mainly 
due to the impossibility for a two-dimensional system to 
be at the same time considered as anisotropic and homo- 
geneous. By contrast, this would be a more reasonable 
assumption in a three dimensional space. Our kinetic 
theory shows too strong mean-field features to be able to 
describe accurately this two-dimensional situation. We 
conclude that the predictions of this isotropic model are 
no longer valid when the domains start growing. 
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tio of the number of growing and shrinking rods n g /n 3 and a 
itself. This Figure illustrates Eq. 1241 



Connections with the individual history of the 
rods 



On the graph showing the individual history of the 
rods (Figure |SJ , one can decompose the time of growth 
T + into a sum of elementary growth intervals t^, and the 
time of shrinking T~ into a sum of elementary shrinking 
intervals t~: 



(17) 



In particular, we expect that the following relations 
between the average elementary time of growth (r + ), 
shrinkage (t~), and the rates p gs and p sg hold for an 
isotropic system: 



Pgs 



(r + >- 



Ps 



(r-y 



(18) 



For an anisotropic system, p gs becomes orientation de- 
pendent, while p sg should not. In a stationary case, the 
following "detailed balance" relation holds: 



^gPgs — T^sPsgi 



(19) 



while in a quasi-stationary state, we expect only a qual- 
itative relation: 

n g ^ Psg 
T^s Pgs 



(20) 



suggesting the relation: 



Psg_ 

Pgs 



(r+) 

(r-y 



(21) 



In any case, from Eq. |7| we have: 

(0 = (v + ~v-)(T-)+v-(T-) 
(T) = <T+) + (T-) . 

The ratio of growing to shrinking time is 



(22) 



(T+) 
(T-) 



av(T) 



1 



w 

v(T) 



(23) 



where the r.h.s. is expected to stay close to the value 
a, as I is smaller than v(T~), except for young rods. 
Because the number of growing time intervals is close 
to the number of shrinking time intervals, there should 
not be a large difference between the ratios (T + ) / (T~) 
and (t + ) /(t~), suggesting, along with Eq. |^ another 
relation: 



(r+> ^ <T+> 
(r-> " (T-> 



(24) 



To test this relation, we performed a set of simula- 
tions with the same rate of injection Qi = 1000, and 



varying a. Both 



"91 



<r+), <t->, (T+) and (T-) 



can 



be independently obtained from the simulation, as sum- 
marized in Figure PHI The ratio n g /n s , T + /T~, and 
(t + ) / (t~) seems to remain remarkably constant as the 
simulation goes on. The agreement is not quantitative, 
but the four quantities of relation|^]show significant cor- 
relations (Figure 



D. Anisotropic and non-stationary solutions of the 
kinetic theory 

There are indications that the stationary, homoge- 
neous situation is not the only possible solution of the 
kinetic theory. We outline in Appendix [Dl the main 
features of a stationary but anisotropic solution, with a 
non trivial dependence of p gs in the orientation Q. This 
solution can explain why the system tends to choose a 
preferential global orientation. However, the predicted 
anisotropy is less than the one that is numerically ob- 
served in our system. 

The kinetic model described above can also describe 
time dependent solutions. However, we did not find 
any simple time-dependent solution compatible with our 
boundary conditions, i.e. a constant injection rate. It 
remains that the existence of an unstable, non station- 
ary and anisotropic solution of the kinetic model cannot 
be ruled out. 



V. CONCLUSION 

We proposed and tested numerically a new paradigm 
for a kinetically constrained growth mechanism of living 
rods. We demonstrated that the collective behavior of 
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the rods were leading to the formation of bundles of well 
oriented rods, in the absence of excluded volume interac- 
tions and chemical gradients. We suggest this mechanism 
as a possible alternative in the formation and the orien- 
tation of microtubules gels. 

We tested the responsive properties of the system to 
some external stress. We found that the alignment was 
much faster in the presence of any small anisotropic bias. 
This result is interesting in the context of the experimen- 
tal microgravity results of Tabony [23, We devel- 
oped a kinetic theory which accounts for the basic scal- 
ing properties of the system: conversion rates, average 
length of the microtubules. The quantitative agreement 
remains poor, due to the strong correlations present in 
this two-dimensional system, in which the emergence of 
large anisotropic domains is incompatible with the ho- 
mogeneity assumption (r independence). 

Because we were able to show that this mechanism also 
induces some alignment in three dimension, we believe 
that our kinetic theory would give a better agreement 
in higher dimensional systems, which is the subject of 
future work. 



APPENDIX A: A KINETIC THEORY FOR 
HOMOGENEOUS DISTRIBUTIONS 

We assume first that the system is homogeneous, and 
that the distribution functions do not depend on r. De- 
noting by S the area of the system, the distribution func- 
tions then reduce to their homogeneous form: 



1 



(Al) 



The functions 



■g,s,t 



are "extensive" functions of the area 



S. We propose a master equation for these functions, 
which accounts for the shrinking and growing behavior 
of the rods, and also account for the possibility of in- 
terconversion between the s and g states. In order to 
cope with a possible global anisotropy of the rods, we let 
the interconversion rates depend on the orientation Q, 
but not on the length: p gs (Q)dt is the fraction of <?-rods 
which switches to a s-rod state during the time inter- 
val dt, while p sg (£l)dt describes the reverse change. The 
corresponding master equation reads: 

c g (7, il, t + dt) = c g (l - vdt, ft, t)+ 

P S g(tt)C S {l, ft, t)dt - Pg S (Cl)Cg(l, O, t)dt] 

c s (l,n,t + dt) =c s (l + v~dt,Sl,t)+ 

P gs (il)Cg(l, fl, t)dt - P sg (fl)Cs(l, O, t)dt. 

The continuous limit dt — » leads to a system of two 
partial differential equations: 



Be 



PsgCs 



Pas 



8Cg_ 

ft* 
oc s oc s 

~~dt ~ aV ~dT = P9sCg ~ PsgCs ' 



(A2) 



The partial derivatives in the left hand sides of Ea. lA2l 
correspond to a drift motion of the rods along the I axis. 
We can associate to this drift the "currents" of grow- 
ing j g (l,Cl,t) — vc g and shrinking j s (l,£l,t) — —avc s 
rods. The sum j g + j s measures the difference between 
the number of rods which have grown bigger than I, and 
the number of rods which have shrunk below I. The dis- 
tribution ct — c g + c s obeys a usual conservation equation 
dct/dt + d(j g + j s )/dl = 0, provided I > 0. 

In particular, as there is no other option for a rod with 
length than growing or disappearing, and our master 
equations must be completed with a boundary condition 
involving j g and the injection rate Qi. 



= vc g Q = 0,fi). 



(A3) 



Meanwhile, the rate of disappearance of the rods is Qd, 
obeying 

Q d (Q) = avc s (l = 0,Sl). (A4) 
At the other extreme, we expect that 

limc ff , s , t (/,fi) = 0. (A5) 

Equations IA2IA3IA5I are the basis of our kinetic theory. 

APPENDIX B: THE STATIONARY CASE 

The above system of equations is simpler if we look 
for a time independent solution, setting the time partial 
derivative to zero. Introducing ct — c g + c s , and = 
c g - c s , we get: 



v ff \ dct , v n , \ dcd n 
v, dc d v dc t 



(Psg -Pgs)Ct ~ (Psg+Pgs)Cd\ 



c(0,ft) = q; lim c(Z,fi) = 0, 

/ — >oo 



with solution 



C t (Z,Q) = qe- 1 ' 1 ^; 
c d (l,n) = 2=±qe- l /m. 



(Bl) 



where, for convenience, we have introduced q — c t (0, f2), 
while the average length in the direction Q is given by: 



1(0) 



ar 



where appears the ratio a = v jv = v /{v + —v ). 



ap gs (tt) -p sg (£Y)' 
Moreover, the distributions c g and c s verify 

Cg (i,n) 



c 8 (i,n) 



: a, 



(B2) 



(B3) 
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and consequently, 



c - _^_ e -'An). 



a+l 



(B4) 



Finally, the injection rate Qi(fl) is related to q(fl) by 



a + l 



(B5) 



where 5 1 is the total area of the system, n t is the total 
number of rods, and (n t /S) is the ratio of two "extensive" 
functions. 



APPENDIX C: THE ISOTROPIC SOLUTION 
AND ITS PREDICTIONS 

The isotropic hypothesis consists in taking p gs and p sg 
angle independent. This simplifies the above kinetic the- 
ory to a point where a self-consistent analytical solution 
becomes available. The f2 dependence of I and l p drops 
out, and we get: 



I = 2 ± 

7T 



which is the basis of a stationary, homogeneous solution 
of the system, expressed in terms of Qi, a, p gs and p sg . 
An example of explicit angular dependence of Qi(£l) is 
the situation described in section Ull Bl and in Figure 01 
In most cases, however, we are interested in an isotropic, 
constant, function Qi, which we consider now. 

To move further, we must estimate the interconversion 
rates p gs and p sg . To estimate p sg , we assume that colli- and equation |B2] leads to the self-consistence relations: 
sions are pairwise, and, as in the usual kinetic theory of 
gases, that there is no correlation between any two col- 
liding rods. Then, p sg does not depend on angles, and is 
inversely proportional to the average waiting time (t~) 
spent in the blocked state. Since all rods shrink at the 
same speed, the waiting time depends only on the dis- 
tance between the contact point and the minus end of 
the restricting rod. Assuming a uniform distribution of 
contact points along the rod, we find that the average 
waiting time associated to a restricting rod with length I 
is l/(2v~), and consequently, 



0" 



+ n t 2l 2v~ 

S 7T J 

3tt S 



2(a+l)n t 



(CI) 



(C2) 



(C3) 



Then, we replace n t by S±1S»_ ) to make a prediction 
for the average length I and the number of rods nt- 



I -\ 1 1 



2v 
1 



(B6) 



I = 



3tt S 



In this equation, we need to know the average length I 
of the rods, irrespective of their orientation. Given the 
solution obtained above, this simply reads: 



I = 



fdn(l(n) 

J d/dfi ct(Z,fi) ~ /dfil(fi) 



2(a + l)n t \n t 

3a Sv f Sv 
4(a + l) 2 Q^ ~ \Q~i 
a+l 



(C4) 

(C5) 
(C6) 



Jdldnict(l,fl) 



(B7) 



The estimate of p gs also comes from the analogy with 
the kinetic theory of gases. We imagine the system from 
the point of view of an observer sitting at the top of a 
growing tip, and estimate the area swept by the mesh of 
all the other rods, moving relatively to the observer at 



speed 



The typical collision time is reached when 



this area becomes comparable to the total area S of the 
system, making the probability of collision of order one. 
The calculation shows that the collision time depends on 
the projected length V of the obstructing rod, and on the 
relative orientation difference Q — fl' between the two 
rods. We can write: 



z p (n) 



Jdfl'dl'l' |sin(ft-fi') | ct(l',Q') 

J dn'di'c t (i',n>) ' 



(B8) 



and the probability p gs (fl) reduces to: 

Pgs (n) = v + r p (n)(^), 



APPENDIX D: THE ANISOTROPIC SOLUTION 

We believe that the system of equations IB2I IB6I IB 71 
IB8IIB9l also admits an anisotropic solution, characterized 
by an explicit Q dependence of p gs (£l) and Z(fi) while p sg 
remains isotropic. To approach this solution, we expand 
1(d) in cosine series: 

!(fi) = l Q + h cos(20) + U cos(4ft) . . . (Dl) 

and approximate sin |fi — f2'| in a similar manner: 

sin = s + s 2 cos(2f2) + s 4 cos(4fi) . . . (D2) 

Possible choices include the Fourier expansion: 

sin |Q| = 1- If (D3) 



P =i 



Ap 2 - 1 ' 



(B9) 

or replacing sin |f2| by sin 2 (0). 
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For instance, by keeping the two first terms in the ex- 
pansion, and writing the self-consistence eauation l_B2l un- 
der the form 1(H) x (ap gs (fl)~p sg ) = av, we finally obtain 
a system of equations for Iq and x — lijl^. 



av + n t lQ 
S 

otv + n t lQ 



so 



S X 



X 

1 + Y 

x" 

1 + T 



s 2 - 



+ s 2 x = 



= 3v~; 
2v~x 



1 



(D4) 



problem. 

Thus, despite its strong mean field features, the kinetic 
model is compatible with the emergence of an anisotropic 
solution, with an explicit angular dependence of the av- 
erage length of the rods. However, this anisotropy is 
bounded, with an average length finite in all directions, 
while the domains observed in the simulations can grow 
without limit. 



We observe that the isotropic solution x = 0, Iq = 
3Sv~ /(s av + n t ) (equivalent to lC3|) coexists along with 
an anisotropic solution x ^ 0, where x solves 



av + n t lf ) 
S 



so 



2tT 



S2 = 



1 T 2 



(D5) 



and Iq is a function of x and the other parameters of the 
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